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Abstract 

In  recent  years,  the  need  for  reliable  modeling  tools  for  lightweight  robotic 
systems  deployed  on  various  terrains  has  spurred  research  efforts  into  develop¬ 
ment  of  vehicle  terrain  interaction  (VTI)  models.  This  paper  presents  an  anal¬ 
ysis  of  rigid  wheels  -  dry  sand  interaction  and  compares  experimental  results 
with  predictions  from  established  terramechanics  theory.  A  novel  experimental 
setup,  based  on  sensing  elements  placed  on  the  wheel  surface,  allows  inference 
of  normal  and  tangential  stress  at  the  wheel-terrain  interface.  A  particle  image 
velocimetry  (PIV)  analysis  is  used  to  study  the  soil  kinematics  under  the  wheel. 
The  analysis  of  stress  profiles  shows  that  stress  patterns  under  lightweight  vehicle 
wheels  conform  reasonably  well  to  established  terramechanics  theory  developed 
for  heavy  vehicles.  For  the  wheel  under  investigation,  the  stress  distribution  had 
minor  variation  along  wheel  width  for  low  slip  conditions.  The  wheel  model 
proposed  by  Wong  and  Reece  was  analyzed  in  light  of  the  stress  and  soil  kine¬ 
matics  measurements  available.  It  was  found  that,  by  appropriately  characteriz¬ 
ing  the  model  coefficients  c  i  and  C2,  and  understanding  the  physical  meaning  of 
the  shear  modulus  kx,  it  is  possible  to  obtain  torque,  drawbar  force,  and  sinkage 
predictions  within  11%  (full  scale  error)  of  experimental  data. 

Keywords:  Wheel  Model,  Shear  Strength,  Stress  Sensor,  Granular  Particle 
Image  Velocimetry,  Terrain  Shearing  Failure,  Wheel  Dynamics,  Off-Road 
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1.  Introduction 

In  recent  years,  the  analysis  of  lightweight  robotic  system  mobility  has  raised 
many  questions  regarding  whether  classical  terramechanics  theory  for  wheeled 
vehicles  is  accurately  predictive  for  reduced  scale  vehicles  [6,  16,  3,  10].  Lack¬ 
ing  a  standardized  classification,  in  this  paper  we  arbitrarily  define  lightweight 
vehicles  as  having  average  ground  pressure  below  20  kPa.  Many  space  rovers 
and  robotic  ground  vehicles  fall  within  this  classification. 

Basing  his  analysis  on  fundamental  concepts  of  soil  mechanics,  Bekker  [2] 
introduced  a  theory  to  predict  mobility  of  wheeled  and  tracked  vehicles  in  off¬ 
road  scenarios.  Bekker  proposed  a  set  of  semi-empirical  equations  to  predict 
different  mobility  aspects,  such  as  compaction  resistance,  traction,  sinkage,  and 
driving  torque.  Bekker  himself  noted  that  terramechanics  theory  “become  less 
accurate  for  wheels  smaller  than  20  inches  [...]  and  for  wheel  loads  below  about 
10  lbs”.  Carrier  [6],  while  studying  the  traffic  ability  of  lunar  micro  rovers,  con¬ 
cluded  that  classical  Bekker  equations  lead  to  an  underestimation  of  small  rover 
tractive  performance.  Richter  et  al.  [16]  investigated  the  performance  of  wheels 
with  diameter  ranging  between  150  mm  and  250  mm  and  vertical  loads  ranging 
from  10  N  up  to  120  N,  and  concluded  that  classical  Bekker  model  needs  correc¬ 
tions  in  order  to  accurately  predict  performance.  Meirion- Griffith  and  Spenko 
[10]  used  small  wheels  as  penetration  plates,  and  noted  that  Bekker’s  pressure- 
sinkage  equation  21  is  affected  by  wheel  curvature.  Griffith  and  Spenko  pro¬ 
posed  a  modified  Bekker  pressure-sinkage  equation  to  account  for  small  wheels’ 
curvature. 

The  theory  for  off  road  rigid  wheel  mobility  evaluation  developed  by  Bekker 
was  further  refined  by  Wong  and  Reece  [23,  24].  Wong  and  Reece  did  not  sim¬ 
ply  apply  correction  factors  to  Bekker  equations,  but  rather  expanded  the  Bekker 
methodology  to  calculate  wheel  performance  through  the  prediction  of  stress 
distributions  at  the  wheel-terrain  interface.  The  model  proposed  by  Wong  and 
Reece  (here  referred  to  as  the  WR  model)  has  the  merit  of  deriving  all  wheel 
performance  metrics  (i.e.,  drawbar  force,  torque,  and  sinkage)  from  the  calcu¬ 
lated  stress  distributions  at  the  interface.  On  the  other  hand,  Bekker’s  original 
approach  was  based  on  a  series  of  ad-hoc  formulations  intended  to  model  each 
single  aspect  of  vehicle  mobility  independently. 

Ishigami  et  al.  [7]  showed  that  a  WR-based  model  could  reasonably  replicate 
single  wheel  experiments  (though  only  positive  slip  was  investigated).  However, 
in  [7]  it  is  not  discussed  how  soil  parameters  were  calculated,  and  therefore  it 
is  reasonable  to  assume  that  some  soil  parameters  were  tuned  to  match  the  ex¬ 
perimental  observations.  Ding  et  al.[3],  noting  a  significant  discrepancy  between 
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measured  and  predicted  sinkage,  proposed  a  modified  WR  model  where  the  sink- 
age  exponent  (see  APPENDIX  A)  is  modified  according  to  slip.  Based  on  the 
authors’  own  experience,  tuning  of  WR  model  parameters  is  inevitably  required 
to  achieve  accurate  model  predictions  across  a  broad  range  of  loading  and  slip 
conditions. 

This  brief  overview  of  the  most  significant  work  on  lightweight  vehicle  mo¬ 
bility  modeling  shows  that  previous  studies  have  either  proposed  modifications 
of  the  Bekker-Wong-Reece  models  (typically  by  introducing  additional  parame¬ 
ters)  or  they  have  arbitrarily  tuned  some  parameters  to  improve  correlation  with 
experimental  data.  In  either  case,  the  inherent  reasons  for  poor  model  perfor¬ 
mance  were  not  investigated. 

To  overcome  these  issues,  in  this  paper  we  describe  a  detailed  analysis  of 
stress  distributions  under  small-sized  rigid  wheels  operating  on  cohesionless  soil 
in  order  to  understand  if,  where,  and  how  WR  models  fail.  (Note  that  the  origi¬ 
nal  Bekker  model  is  not  discussed).  A  custom  force  sensing  array  located  at  the 
wheel-terrain  interface  is  used  to  measure  stresses  at  the  wheel  interface.  The 
force  sensors  are  strain  gage-based  flexural  elements  with  interchangeable  inter¬ 
face  surfaces  that  are  designed  for  integration  with  wheels  or  other  running  gear. 
The  sensors  allow  explicit  measurement  of  normal  and  shear  forces  (and,  there¬ 
fore,  estimation  of  normal  and  shear  stresses)  at  numerous  discrete  points  along 
the  wheel-soil  interface.  Similar  experimental  methodologies  were  employed  by 
Hegedus  [4],  Sela  [17],  Onafeko  and  Reece  [15],  Krick  [9],  and  Shamay  [20]. 
The  key  difference  is  that  in  [4,  17,  15,  9,  20]  the  average  wheel  ground  pressure 
was  approximately  100  kPa,  while  in  this  paper  the  wheel  average  ground  pres¬ 
sure  is  on  the  order  of  10  kPa.  (The  average  ground  pressure  is  evaluated  as  the 
nominal  wheel  load  distributed  over  a  flat  wheel  section  spanning  30  degrees). 
Oida  et  al.  [14]  instrumented  a  flexible  tire  with  a  sensor,  based  on  Krick’s  de¬ 
sign  [9],  and  they  measured  normal,  tangential,  and  lateral  stress  at  the  tire-sand 
interface  (however  the  vertical  load  and  tire  dimensions  are  unknown).  Nagatani 
et  al.  [13]  have  used  stock  button-type  force  transducers  to  measure  normal  stress 
at  wheel-terrain  interface.  Although  the  average  ground  pressure  was  compara¬ 
ble  to  what  is  studied  here,  the  setup  proposed  in  [13]  was  only  able  to  measure 
normal  load. 

Another  experimental  methodology  employed  in  this  work  relies  on  imaging 
of  the  wheel-soil  interface  and  the  use  of  particle  image  velocimetry  (PIV)  to 
measure  micro-scale  terrain  displacement.  This  methodology,  although  confined 
to  a  plane  strain  case,  allows  measurement  of  the  soil  displacement  field  under 
the  wheel.  Though,  this  method  does  not  explicitly  permit  calculation  of  the 
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Table  1:  Wong  and  Reece  wheel  model  terrain  parameters  and  coeffi¬ 

cients. 


Symbol 

Units 

Description 

n 

n/a 

sinkage  exponent 

kc 

kN/m'!+1 

pressure-sinkage  coefficient 

k(j) 

kN/m'l+2 

pressure-sinkage  coefficient 

c 

Pa 

cohesion 

§ 

deg 

angle  of  internal  friction 

K 

m 

shear  deformation  modulus 

Cl  ,C2 

n/a 

coefficients  for  determining  the  relative  position 
of  maximum  radial  stress 

0, 

deg 

exit  angle 

velocities  of  individual  soil  particles,  it  does  allow  estimation  of  a  regularly- 
spaced  velocity  field  in  the  soil.  While  such  visualization  techniques  have  been 
widely  employed  in  the  field  of  experimental  fluid  mechanics,  their  application 
to  the  study  of  soils  is  a  relatively  new  development  [21,  11,  12]. 

Measurements  of  stress  distributions  and  the  soil  velocity  field  are  comple¬ 
mented  by  an  in-depth  comparison  with  WR  model  predictions.  The  model  relies 
on  a  set  of  6  terrain  parameters  and  3  wheel-terrain  interaction  coefficients,  pre¬ 
sented  in  Table  1.  This  work  identifies  the  shear  modulus,  kx,  and  the  coefficients 
for  determining  the  relative  position  of  the  maximum  radial  stress,  c\  and  c2,  as 
the  principal  factors  that  often  lead  to  poor  performance  of  WR  model  predic¬ 
tions. 

Here,  we  have  confined  our  study  to  wheel  operation  on  dry  sand.  The  sand 
utilized  in  this  paper  has  been  fully  characterized  via  a  series  of  direct  shear  tests 
(ASTM  D3080)  and  penetration  tests.  Direct  shear  tests  were  performed  to  esti¬ 
mate  shearing  parameters  such  as  cohesion,  c,  angle  of  internal  friction,  (j),  and 
shear  modulus  kx.  Penetration  tests,  although  not  standard  tests,  were  performed 
to  evaluate  the  “Bekker  parameters”  n,  kc,  and  k^,  which  are  necessary  for  char¬ 
acterization  of  the  pressure-sinkage  behavior  of  the  soil.  The  key  questions  that 
this  paper  addresses  are: 

Ql  Are  the  stress  distribution  that  form  under  lightweight  wheels  similar  in  na¬ 
ture  to  those  that  form  under  heavy  weight  vehicles? 

Q2  Is  the  WR  wheel  model  capable  of  accurately  modeling  lightweight  vehicle 
mobility? 
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The  paper  is  organized  as  follows:  Experimental  Setup  presents  the  terrame- 
chanics  rig,  the  custom  sensing  array,  and  the  PIV  setup  utilized  for  this  work. 
Experimental  Data  Collection  and  Discussion  section  is  articulated  in  6  subsec¬ 
tions:  first  we  present  an  overview  of  measured  stress  distributions  under  the 
wheel;  then  we  compare  the  output  produced  by  the  custom  sensing  array  with 
the  readings  obtained  by  other  sensors  installed  on  the  wheel/rig;  subsequently 
we  investigate  coefficients  c\  and  C2  and  how  these  affect  wheel  slip-sinkage  be¬ 
havior;  we  then  proceed  to  analyze  PIV  data  in  order  to  understand  soil  shear 
displacement  modelling  at  the  wheel-soil  interface;  next  we  use  this  analysis  to 
understand  how  the  shear  modulus  kx  affects  shear  stress;  finally  we  present  a 
modified  WR  model. 

2.  Experimental  Setup 

2.1.  Single  Wheel  Test  Rig 

The  Robotic  Mobility  Group  at  MIT  has  designed  and  fabricated  a  multipur¬ 
pose  terramechanics  rig  based  on  the  standard  design  described  by  Iagnemma 
[5].  The  testbed  is  pictured  in  Figure  1(a)  and  it  is  composed  of  a  Lexan  soil 
bin  surrounded  by  an  aluminum  frame  where  all  the  moving  parts,  actuators  and 
sensors  are  attached.  A  carriage  slides  on  two  low-friction  rails  to  allow  longi¬ 
tudinal  translation  while  the  wheel  or  track,  attached  to  the  carriage,  is  able  to 
rotate  at  a  desired  angular  velocity.  The  wheel  mount  is  also  able  to  freely  trans¬ 
late  in  the  vertical  direction.  This  typical  setup  allows  control  of  slip  and  vertical 
load  by  modifying  the  translational  velocity  of  the  carriage,  angular  velocity 
of  the  wheel,  and  applied  load.  Horizontal  carriage  displacement  is  controlled 
through  a  toothed  belt  actuated  by  a  90  W  Maxon  DC  motor,  while  the  wheel  is 
directly  driven  by  a  200  W  Maxon  DC  motor.  The  motors  are  controlled  through 
two  identical  Maxon  ADS  50/10  4-Q-DC  servoamplifiers.  The  carriage  hori¬ 
zontal  displacement  is  monitored  with  a  Micro  Epsilon  WPS-1250-MK46  draw 
wire  encoder  while  wheel  vertical  displacement  (i.e.,  sinkage)  is  measured  with 
a  Turck  A50  draw  wire  encoder. 

A  6-axis  force  torque  ATI  Omega  85  transducer  is  mounted  between  the 
wheel  mount  and  the  carriage  in  order  to  measure  vertical  load  and  traction  gen¬ 
erated  by  the  wheel.  Finally,  a  flange-to-flange  reaction  torque  sensor  from  Futek 
(TFF500)  is  used  to  measure  driving  torque  applied  to  the  wheel.  Control  and 
measurement  signals  are  handled  by  a  NI  PCIe-6363  card  through  Lab  view  soft¬ 
ware.  The  rig  is  capable  of  approximately  1  meter  of  horizontal  displacement  at 
a  maximum  velocity  of  approximately  120  mm/s  with  a  maximal  wheel  angular 
velocity  of  approximately  40  deg/s.  The  bin  width  is  0.6  meters  while  the  soil 
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Horizontal 


(b) 


Figure  1:  (a)  The  terramechanics  rig  at  MIT  (b)  equipped  with  imager  and  lighting  for  PIV 
analysis. 

depth  is  0.16  meters.  Considering  the  wheel  sizes  and  vertical  loads  under  study, 
these  physical  dimensions  are  sufficient  for  eliminating  boundary  effects. 

Moreover,  the  same  testbed,  with  some  adaptations,  can  be  used  to  perform 
soil  penetration  tests.  For  the  experiments  described  in  this  paper,  the  Mojave 
Martian  Simulant  (MMS)  was  employed  as  a  test  medium  [1],  MMS  is  a  mixture 
of  finely  crushed  and  sorted  granular  basalt  intended  to  mimic,  both  at  chemical 
and  mechanical  levels,  the  Mars  soil  characteristics.  MMS  particle  size  distribu¬ 
tion  spans  from  micron  to  millimeter  scale,  with  80%  of  particles  above  the  10 
micron  threshold.  Soil  properties  were  measured  through  a  series  of  plate  pene¬ 
tration  tests  and  direct  shear  tests:  nominal  MMS  soil  parameters  are  presented 
in  Table  2.  It  should  be  noted  that  the  shear  modulus  is  very  small.  Typical  liter¬ 
ature  values  range  between  0.01  and  0.03  m,  however,  as  was  presented  in  [18], 
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Table  2:  Mojave  Martian  Simulant  (MMS)  properties  measured  through  a  series  of  plate  penetra¬ 
tion  tests  and  direct  shear  tests.  It  should  be  noted  that  the  shear  modulus  is  very  small.  Typical 
literature  values  range  between  0.01  and  0.03  m,  however,  as  was  presented  in  [18],  correct  calcu¬ 
lation  of  shear  modulus  leads  to  a  significantly  smaller  value  of  kx. 


Symbol 

Value 

Units 

n 

1.4 

n/a 

kc 

846 

kN/m"+1 

k(j> 

6708 

kN/m"+2 

c 

600 

Pa 

0 

35 

deg 

0.0006 

m 

correct  calculation  of  shear  modulus  leads  to  a  significantly  smaller  value  of  kx. 
2.2.  Force  Transducers 

The  sensing  elements  used  to  measure  forces  at  the  wheel-terrain  interface 
are  custom-made  flexural  devices  equipped  with  strain  gages;  similar  designs 
have  been  used  in  [9,  20].  Each  sensor  is  composed  of  an  L-shaped  beam 
equipped  with  350  ohm  strain  gages  placed  on  the  two  arms  of  the  fixture  as 
shown  in  Figure  2(a)  .  Two  pairs  of  strain  gages  are  connected  to  each  surface 
of  the  beam  element  in  a  full-bridge  configuration,  allowing  maximal  bending 
response  and  naturally  rejecting  axial  loading.  Considering  the  symmetry  of  the 
problem,  only  half  of  the  wheel  width  was  equipped  with  sensors,  as  shown  in 
Figure  2(b). 

The  sensor  was  designed  to  maximize  strain  at  points  g\  and  g2  while  keep¬ 
ing  the  sensor  head  displacements  below  0.2  mm.  Considering  0-100  kPa  as  a 
pressure  range,  the  sensor  was  designed  to  withstand  loads  up  to  18  N,  with  a 
sensitivity  of  approximately  0.01  N.  If  F\  and  /3  are  the  tangential  and  normal 
loads  acting  at  the  tip  of  the  sensor,  gi  and  g2  are  the  readings  from  the  strain 
gages,  d\  and  di  are  the  distances  of  the  gages  from  the  tip,  and  £  is  the  normal 
load  offset;  the  following  system  of  linear  equations  holds: 


g  t  =  F\d\  +F2£ 
g2  —  F\d\  +  F2J2  +  F2& 

Solving  for  F\  and  Fr. 


(1) 

(2) 


F2 


82- gi 
d2 


7 


(3) 
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(a)  (b) 


Figure  2:  (a)  Schematic  of  the  custom  force  sensor  for  interfacial  stress  measurement.  Sensor 
body  is  12  mm  wide  while  sensor  head  has  an  area  of  10  mm  x  18  mm.  (b)  Five  sensors  are 
distributed  from  the  wheel  median  axis  to  the  wheel  edge.  Sensors  are  rigidly  connected  to  the 
wheel  hub. 


gl  F2E 
d\  d\ 


(4) 


Therefore,  the  tangential  load  F\  is  not  independent  of  Fi  in  the  presence 
of  an  asymmetric  pressure  distribution  over  the  sensor’s  head.  (If  the  pressure 
distribution  is  uniform,  the  offset  £  would  be  zero).  Assuming  that  the  normal 
pressure  distribution  is  linear  along  the  sensor  head,  the  offset  £  is  equal  to  ~  1.6 
mm,  and  therefore  F\  can  be  corrected  once  Fi  is  calculated.  It  should  be  noted 
that  the  offset  error  is  minimized  by  design,  because  £/d\  ~=  0.03,  while  the 
normal  stress  distribution  is  non-uniform  only  in  the  proximity  of  the  terrain  en¬ 
try  and  exit  regions  (where,  on  the  other  hand,  normal  stress  is  small  in  absolute 
value). 

Calibration  of  the  flexure  elements  was  performed  by  loading  the  sensors  in 
the  tangential  and  normal  direction  with  test  weights  of  100  g,  200  g,  and  500 
g.  The  results  (presented  in  Table  3)  show  that  the  transducer  error  is  always 
below  2%.  Calibration  factors  were  obtained  by  minimizing  the  least  square 
error  among  all  the  calibration  tests  performed.  When  assembled  on  the  wheel, 
the  sensing  array  showed  to  be  extremely  sensitive  to  transducer  positioning: 
minimal  misalignment  of  transducer  heads  would  produce  skewed  readings  due 
to  the  uneven  contact  geometry.  This  issue  was  controlled  by  repeatedly  testing 
the  wheel  on  a  hard,  flat  surface,  and  adjusting  individual  sensor  positions  while 
verifying  a  uniform  output  across  all  sensors. 
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Table  3:  Force  sensing  array  calibration.  Accuracy  is  within  2%. 


Fi 

Sensor  I 

Sensor  II 

Sensor  III 

Sensor  IV 

Sensor  V 

Nominal 

Reading 

Error 

Reading 

Error 

Reading 

Error 

Reading 

Error 

Reading 

Error 

Load  [g] 

[g] 

[%] 

[g] 

[%] 

[g] 

[%] 

[g] 

[%] 

[g] 

[%] 

100 

100.1 

0.14 

101.0 

0.99 

99.7 

-0.29 

100.1 

0.09 

100.1 

0.14 

200 

201.7 

0.84 

200.8 

0.42 

200.6 

0.29 

200.2 

0.12 

196.8 

-1.58 

500 

498.9 

-0.22 

498.9 

-0.21 

499.7 

-0.07 

499.5 

-0.11 

500.7 

0.13 

Fi 

Sensor  I 

Sensor  II 

Sensor  III 

Sensor  IV 

Sensor  V 

Nominal 

Reading 

Error 

Reading 

Error 

Reading 

Error 

Reading 

Error 

Reading 

Error 

Load  [g] 

[g] 

[%] 

[g] 

[%] 

[g] 

[%] 

[g] 

[%] 

[g] 

[%] 

100 

98.2 

-1.76 

100.8 

0.82 

100.0 

0.00 

98.5 

-1.49 

101.3 

1.28 

200 

199.1 

-0.45 

199.9 

-0.04 

200.6 

0.30 

199.2 

-0.42 

201.2 

0.62 

500 

499.5 

-0.10 

499.8 

-0.05 

499.8 

-0.05 

500.4 

0.07 

500.5 

0.10 

2.3.  Particle  Image  Velocimetry 

For  PIV  experiments,  the  Lexan  soil  bin  was  fitted  with  a  0.0254  m  thick  tem¬ 
pered  glass  wall  while  the  running  gear  was  operated  flush  against  this  surface 
(see  Figure  1(b)).  It  should  be  noted  that  stress  measurements  and  particle  image 
velocimetry  experiments  were  conducted  independently  (i.e.,  a  wheel  with  the 
exact  same  dimensions,  but  not  equipped  with  interfacial  force  transducers,  was 
used  for  PIV  experiments). 

Image  sets  for  the  PIV  measurements  were  captured  with  a  Phantom  7  high¬ 
speed  camera.  The  Phantom  7  is  able  to  record  grayscale  images  at  a  maximum 
resolution  of  800x600  pixels  at  a  maximum  frame  rate  of  6688  fps.  The  camera 
was  placed  perpendicular  to  the  front  glass  wall  (see  Figure  1(b))  at  a  distance 
of  0.52  m,  while  its  focal  length  was  set  to  77  mm  (a  zoom  lens  was  used)  re¬ 
sulting  in  an  image  capture  region  of  approximately  0.15  x  0.11  m.  It  should 
be  noted  that  determination  of  image  capture  region  size  is  largely  dictated  by 
the  particular  experimental  conditions.  Here,  the  image  capture  region  was  cho¬ 
sen  to  conservatively  bound  the  region  of  soil  that  would  undergo  motion  when 
subjected  to  wheel  passage  on  the  soil  surface.  As  noted  previously,  the  MMS 
particle  size  distribution  spans  from  the  micron  level  to  mm  level  with  80  %  of 
particles  above  the  10  micron  threshold.  For  the  imager  configuration  described 
below,  the  average  particle  density  resulted  in  approximately  0.044  pixels  per 
particle. 

Since  soil-glass  friction  could  not  be  accurately  controlled,  it  should  be  noted 
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that  calculated  velocities  and  strain  are  influenced  by  soil-glass  friction.  When 
setting  up  a  similar  experiment,  Wong  [22]  proposed  to  use  a  wheel  with  half  the 
nominal  width,  and  to  reduce  the  normal  load  to  half  as  well.  The  rationale  was 
that,  assuming  negligible  frictional  drag,  the  glass  becomes  a  plane  of  symmetry, 
and  therefore  the  soil  motion  is  analogous  to  the  motion  under  the  median  axis 
of  a  full  sized  wheel. 

In  this  paper,  we  did  not  follow  this  approach  for  two  reasons:  vertical  load 
and  wheel  size  do  not  scale  linearly;  and  even  assuming  linear  scaling  this  would 
require  wheels  of  different  width  for  every  vertical  load  to  be  tested.  Therefore, 
assuming  limited  soil-glass  friction  and  minor  stress  non-uniformity  along  the 
wheel  width  [19],  the  soil  motion  at  the  glass  interface  remains  representative  of 
the  soil  motion  under  the  wheel. 

3.  Experimental  Data  Collection  and  Discussion 

Tests  were  conducted  with  a  smooth  aluminium  wheel  of  0.13  m  radius  and 
0.16  m  width.  The  wheel  was  coated  with  a  layer  of  soil  in  order  to  ensure 
an  adequate  level  of  friction  between  the  wheel  surface  and  the  terrain.  This 
wheel  has  approximately  the  same  dimension  of  a  NASA  Mars  exploration  rover 
(MER),  a  successful  lightweight  robotic  system.  Three  vertical  loads  have  been 
investigated,  70  N,  100  N  (nominal  for  MER  rovers),  and  150  N,  while  9  slip 
levels  were  selected:  -70%,  -50%,  -30%,  -10%,  0%,  +10%,  +30%,  +50%,  +70%. 
During  experiments,  the  angular  velocity  of  the  wheel  was  held  constant  (17 
deg/s),  while  longitudinal  velocity  was  varied,  for  each  slip  level,  according  to 
the  following  equation: 


where  v  is  the  wheel  longitudinal  velocity,  r  is  wheel  radius,  and  a)  is  wheel 
angular  velocity. 

Note  that  the  same  definition  of  slip  was  used  for  positive  and  negative  slip 
tests.  Each  test  was  repeated  at  least  15  times,  and  the  data  presented  here  was 
obtained  as  the  average  of  all  the  trials. 

Images  for  PIV  analysis  were  collected  while  running  the  wheel  flush  against 
the  2.54  cm  thick  glass  wall.  For  these  experiments,  a  wheel  was  used  with  the 
exact  same  dimension  but  without  a  slot  opening  for  the  force  sensors. 

The  ability  to  measure  normal  and  tangential  forces  at  the  wheel-terrain  in¬ 
terface,  and  the  soil  kinematics  under  the  wheel,  presents  an  opportunity  to  in¬ 
vestigate  the  validity  of  the  WR  model  when  applied  to  lightweight  vehicles. 
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This  model  is,  in  fact,  based  on  the  estimation  of  stress  profiles  along  the  wheel 
contact  patch,  and  therefore  it  is  intimately  dependent  on  the  ability  to  properly 
predict  normal  and  tangential  stresses  along  the  wheel-terrain  interface.  A  de¬ 
tailed  description  of  the  model  is  given  in  APPENDIX  A;  here  only  the  equations 
relevant  to  this  analysis  will  be  discussed. 

This  section  is  organized  as  follows.  A  qualitative  overview  of  measured 
stress  data  is  presented  first.  Then,  the  force  sensor  reading  is  compared  to  that 
from  an  ATI  force  sensor  and  Futek  torque  sensor;  this  allows  us  to  understand 
if  the  WR  approach  (i.e.,  calculate  stress  profiles  and  derive  drawbar,  torque,  and 
sinkage  from  there)  is  valid  for  lightweight  vehicles.  Next,  it  is  shown  how  the 
model  parameters  c \  and  C2  are  calculated  from  stress  measurements.  Finally, 
an  analysis  of  the  shear  modulus  kx  is  performed  and  a  modified  WR  model  is 
presented. 

3.1.  Wheel-Terrain  Interfacial  Stress 

Stress  readings  collected  at  various  vertical  loads  were  observed  to  follow 
similar  trends,  and  therefore  only  data  for  Fz  -  100  N  is  presented  in  Figure  3. 
These  data  illustrate  how  normal  and  tangential  stress  distributions  vary  along 
wheel  width  under  low  normal  load.  For  low  slip,  the  variation  is  minimal,  while 
for  high  slip,  stresses  at  the  wheel  edges  (for  high  positive  slip)  or  at  the  middle 
(for  high  negative  slip)  increase.  Results  for  Fz  =  100  N  are  summarized  in  Table 
4  where  the  peak  normal  stress  for  each  slip  level  is  reported. 

We  hypothesize  that  stress  variation  at  high  slip  is  caused  by  the  soil  flow 
transport  under  the  wheel:  for  high  positive  slip,  the  wheel  displaces  soil  pref¬ 
erentially  along  the  median  axis.  This  is  because  the  soil  at  either  wheel  edge 
remains  stationary,  and  therefore  material  is  preferentially  transported  through 
the  middle.  (The  hypothesized  soil  motion  profile  is  qualitatively  similar  to  the 
fluid  velocity  profile  present  in  open  channel  flow,  where  fluid  velocity  has  a 
peak  value  at  the  channel  centerline  and  zero  velocity  at  the  channel  edges).  As 
a  result,  the  wheel  load  is  preferentially  supported  by  the  wheel  edges,  where 
relatively  less  soil  is  displaced. 

For  negative  slip,  wheel  motion  tends  to  bulldoze  material  at  the  wheel’s 
fore.  For  the  same  reasons  expressed  for  the  high  positive  slip  case,  material  is 
transported  preferentially  along  the  centerline,  where  it  tends  to  accumulate.  As 
a  result,  the  wheel  travels  on  top  of  this  small  mass  of  material  which,  being  at 
higher  elevation  along  the  centerline,  increases  stress  along  the  wheel  median 
axis. 

It  should  be  noted  that  Krick  [9]  reports  opposite  behavior  for  positive  slip 
(negative  slip  is  not  tested):  stress  at  the  wheel  edges  decreases  for  high  slip. 
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(C) 


(d) 


Figure  3:  Normal  and  tangential  stress  for  Fz  =  100  N.  (a)  -30%  slip,  (b)  +30%  slip,  (c)  -70% 
slip,  (d)  +70%  slip.  These  plots  were  obtained  by  averaging  over  all  the  tests  (15  repetitions 
were  conducted  for  each  slip  and  load  level).  The  boxplot  represents  the  standard  deviation  for 
the  data  point. 
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Table  4:  Peak  normal  stress  for  Fz  =  100  N.  Readings  from  sensors  II,  III,  IV,  and  V  are  shown 
in  absolute  value  and  as  percent  variation  with  respect  to  sensor  I.  Large  variations  are  observed 
for  sensors  IV  and  V  at  high  slip. 


i 

o1 

max 

[kPa] 

®max 

[kPa]  [%] 

Gmax 

[kPa]  [%] 

o 

[kPa] 

IV 

max 

[%] 

max 

[kPa]  [%] 

-0.7 

22.8 

20.5 

-10.1 

20.7 

-9.0 

19.8 

-12.9 

15.9 

-30.3 

-0.5 

24.6 

23.1 

-6.3 

22.6 

-8.1 

21.1 

-14.2 

16.7 

-32.3 

-0.3 

28.5 

27.2 

-4.6 

26.2 

-8.0 

25.1 

-12.1 

20.4 

-28.5 

-0.1 

27.3 

26.5 

-2.8 

25.8 

-5.4 

26.3 

-3.5 

24.7 

-9.6 

0 

26.4 

27.1 

2.7 

24.7 

-6.6 

25.0 

-5.4 

25.3 

-4.1 

0.1 

24.5 

24.5 

-0.2 

23.4 

-4.7 

22.6 

-7.8 

19.8 

-19.0 

0.3 

18.4 

17.9 

-3.1 

17.9 

-3.1 

20.5 

11.2 

19.9 

8.1 

0.5 

17.3 

17.3 

0.3 

17.2 

-0.4 

19.7 

13.9 

19.0 

10.3 

0.7 

14.3 

14.2 

-1.0 

14.6 

2.2 

19.7 

37.4 

18.9 

32.3 

Excluding  systematic  experimental  error,  the  reason  could  be  explained  by  the 
different  aspect  ratio  of  the  wheels:  in  this  paper  w/r  >  1  while  in  [9]  w/r  <  1. 
(Unfortunately  [9]  does  not  report  the  wheel  dimension;  however,  from  pho¬ 
tographic  evidences  it  is  clear  that  wheel  width  is  smaller  than  wheel  radius). 
Interestingly,  Shamay  [20],  whose  experiments  are  similar  to  [9],  does  not  report 
higher  stress  at  the  edge  of  the  wheel.  On  the  contrary,  in  Shamay’s  experiments, 
the  stress  at  the  median  axis  of  the  wheel  is  always  dominant.  Hegedus  [4]  only 
measured  normal  stress  at  three  locations  (median  axes  and  both  wheel  edges) 
and  for  positive  slip  and  he  noted  that  stress  at  the  wheel  median  axis  was  always 
dominant.  It  should  be  mentioned  that  peak  normal  stress  in  the  aforementioned 
studies  was  always  above  50  kPa. 

On  the  other  hand,  the  wheel  dimensions  and  loads  used  in  Nagatani  et  al. 
[13]  experiments  are  closer  to  the  ones  utilized  in  this  work,  and  they  suggest 
that  normal  pressure  is  dominant  at  the  wheel  median  axis.  However,  it  should 
be  noted  that  in  [13]  only  one  slip  level  is  presented,  and  shear  stress  is  not 
measured  at  all.  Moreover,  a  load  button  type  load  cell  is  used  with  only  four 
sensing  elements  across  the  entire  wheel  width. 

In  experiments  by  Sela  [17],  Onafeko  and  Reece  [15],  and  Oida  et  al.  [14] 
the  sensing  device  was  not  sensitive  to  stress  variation  across  the  wheel  width, 
therefore  it  is  not  possible  to  draw  any  conclusion,  from  these  studies. 
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3.2.  Comparison  of  Stress  Readings  with  Conventional  Force/Torque  Sensor 

In  order  to  evaluate  the  accuracy  of  the  sensing  array,  and  to  verify  the  WR 
methodology,  it  is  interesting  to  compare  readings  from  the  custom  sensing  array 
with  ones  derived  from  other  sensors. 

The  ATI  force-torque  sensor  continuously  measures  vertical  load  Fz  and  trac¬ 
tive  force  Fx  acting  on  the  wheel.  (The  load  cell  measures  forces  and  torques  in 
all  directions,  but  only  the  aforementioned  quantities  are  relevant  to  this  analy¬ 
sis).  Therefore  it  is  possible  to  compare  the  readings  obtained  by  the  ATI  load 
cell  with  the  stress  measurements  from  the  sensing  array  as  follows: 


FATI 


Of 


wr  I  TSAsin  0d6  +  wr  I  aMcos  BdO 
J  0.  J  Of 


rO 


7 


SA, 


(6) 


r  Of  rOf 

F^TImwr  iSA  cos  0d6  —  wr  /  cT^sinfldfl  (7) 

J  0r  J  0r 

where  FAT1  and  FATI  are  the  vertical  load  and  drawbar  force  as  measured  by  the 
ATI  load  cell,  while  zSA  and  oSA  are  the  tangential  and  normal  stress  as  measured 
by  the  custom  sensing  array.  The  angles  Of  and  0,  are  derived  from  the  sensing 
array  readings,  while  w  and  r  are  respectively  the  wheel  width  and  radius. 

Similarly,  the  Futek  torque  sensor  continuously  measures  the  torque  applied 
to  the  wheel  by  the  motor.  Thus,  it  is  possible  to  compare  the  sensing  array 
readings  as  follows: 

MFutek  fa  wr2  f  f  r SAdQ  (8) 

hr 

where  MFutek  is  the  torque  applied  to  the  wheel  axle,  as  measured  by  the  Futek 
torque  sensor. 

Although  the  test  rig  is  equipped  with  a  draw  wire  encoder  that  measures  the 
vertical  position  of  the  wheel  assembly,  measurements  from  this  sensor  tend  to 
be  very  sensitive  to  soil  preparation  (i.e.,  uneven  terrain  surface).  For  this  reason 
sinkage  was  measured  adopting  a  different  methodology.  Using  a  caliper,  the 
depth  of  the  rut  left  by  the  wheel  was  measured  after  each  test  was  completed. 
The  resulting  wheel  sinkage  can  be  compared  with  the  sensing  array  readings 
through  the  following  equation: 


z 


caliper 


r{  1  —  COS0^A) 


(9) 


where  zcaliper  is  the  sinkage  manually  measured  with  the  caliper  and  0AA  is  the 
entry  angle  derived  from  the  sensing  array  readings. 
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Figure  4:  Readings  from  the  custom  sensing  array  are  compared  against  readings  obtained  from 
the  other  sensors  installed  on  the  wheel.  Comparison  of  (a)  vertical  load,  (b)  torque,  (c)  sinkage, 
and  (d)  tractive  force  are  presented.  Data  is  collected  at  9  discrete  points  (see  Table  5)  but  here 
is  presented  as  continuous  lines  to  improve  visualization. 


Data  presented  in  Figure  4  show  that  the  sensing  array  is  in  good  agreement 
with  data  measured  through  other  sources.  Because  of  an  unrecoverable  sensor 
failure,  drawbar  pull  data  (i.e.,  Fx)  for  Fz  =  1 50  N  is  not  available.  Torque  reading 
was  not  available  for  Fz  -  150  N  and  slip  levels  above  30%  because  of  sensor 
saturation. 

Table  5  shows  the  variation  between  the  stress  array  readings  and  readings 
from  the  other  instruments.  Two  error  measures  are  presented:  the  absolute  error 
A  and  the  full  scale  percent  error  £j. 

A  =  XSA  -  X other  (10) 


£/  = 


15 


A 
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Table  5:  Error  between  the  sensing  array  readings  and  the  other  sensors  mounted  on  the  wheel. 
The  full  scale  error  Ef  shows  that  the  sensing  array  readings  and  other  sensors  readings  are  in 
good  agreement.  The  largest  deviations  are  observed  for  drawbar  force  Fx  readings. 
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where  X$a  is  the  measure  (either  vertical  load  Fz,  traction  Fx,  torque  M,  or  sink- 
age  z)  obtained  through  the  custom  sensing  array,  Xot/7er  is  the  same  quantity 
obtained  with  the  ATI,  Futek,  or  caliper  sensor,  and  Xf  is  the  full  scale  value  for 
the  quantity  under  investigation.  The  full  scale  values  were  selected  assuming 
maximum,  realistic,  operation  conditions:  Fz  =  150  N,  Fx  =  50  N,  M  -  15  Nm,  z 
=  30  mm. 

Overall,  Fx  is  the  measure  that  differs  most  on  average,  while  more  than  half 
of  the  data  points  are  below  5%  full  scale  percent  error  £/. 

3.3.  Calculation  of  Coefficients  for  Determining  the  Relative  Position  of  Maxi¬ 
mum  Radial  Stress 

The  WR  model  [23,  24]  depends  on  6  terrain  parameters  (see  Table  1)  and 
three  coefficients  (Gr,  c\ ,  cf)-  While  the  exit  angle  is  usually  approximated  within 
the  range  of  -10  to  0  degrees,  there  is  virtually  no  knowledge  of  how  the  coeffi¬ 
cients  ci  and  C2  can  be  estimated.  The  WR  model  in  fact  predicts  that  the  normal 
stress  will  reach  a  maximal  value  for  an  angle  9m  that  can  be  calculated  as: 

0m  =  (ci  +c2i)0f  (12) 

This  linear  relationship  was  introduced  by  Wong  and  Reece  [23,  24]  based  on 
the  observations  of  Onafeko  and  Reece  [15].  Unfortunately,  to  date  its  accuracy 
has  not  been  verified  for  lightweight  vehicles.  Also,  it  is  not  known  if  these 
parameters  are  independent  of  wheel  vertical  load.  The  custom  sensing  array 
enables  us  to  investigate  this  matter  further,  since  we  can  directly  measure  the 
angular  location  of  maximum  normal  stress. 

Figure  5(a)  presents  (0y,  9m.  Of)  as  functions  of  slip,  for  all  tested  load  levels. 
In  general,  9m  grows  approximately  linearly  as  positive  slip  increases.  On  the 
other  hand,  for  negative  slip,  9in  stays  relatively  constant. 


Table  6:  Coefficients  for  determining  the  relative  position  of  maximum  normal  stress.  The  table 
includes  the  values  calculated  by  this  study  (second  and  third  column)  together  with  the  data 
available  in  the  literature.  Oida  et  al.  [14]  presented  a  set  of  coefficients,  however  these  are 
omitted  here  because  they  were  obtained  with  a  flexible  tire. 


Coefficient  Positive  Negative  Positive  Slip  Positive  Slip  Positive  Slip  Positive  Slip 
Slip  MIT  Slip  MIT  Onafeko  I  [23]  Onafeko  II  [23]  Hegedus  [23]  Sela  [23] 
d  0.38  0.35  0.43  0.18  0.285  0.38 

c2  0.44  0.11  0.32  0.32  0.32  0.41 
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Figure  5:  (a)  Contact  angles  for  Fz  =  70  N,  Fz  =  100  N,  Fz  =  150  N,  as  function  of  wheel  slip,  (b) 
Slip  sinkage  characteristics  obtained  using  Wong  and  Reece  model  and  parameters  presented  in 
Table  6.  This  shows  how  coefficients  c \  and  c 2  influence  slip-sinkage  behavior. 
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Another  interesting  observation  is  that  the  exit  angle,  Q,  follows  a  similar 
trend:  it  stays  approximately  constant  for  negative  slip  ratios  and  grows  for  pos¬ 
itive  slip  ratios. 

Constants  c\  and  C2,  for  both  positive  and  negative  slip,  were  obtained  by 
linear  regression  and  are  presented  in  Table  6,  together  with  the  data  available  in 
the  literature.  Since  it  is  believed  that  coefficients  depend  on  wheel-terrain  con¬ 
figuration,  it  is  difficult  to  directly  compare  data  because  wheel  size  and  terrain 
type  utilized  in  this  work  is  different  from  [4,  17,  15]. 

Nonetheless,  it  should  be  noted  that  c\  and  C2  control  slip-sinkage  behavior, 
and  therefore,  when  they  are  not  explicitly  available,  they  can  be  used  as  tuning 
parameters  in  order  to  obtain  desired  slip-sinkage  characteristics.  Ding  et  al.  [3] 
suggested  modification  of  the  sinkage  exponent  n  to  improve  slip-sinkage  model 
response.  However,  it  is  usually  not  possible  to  directly  measure  the  coefficients 
ci  and  C2  (as  opposed  to  the  sinkage  exponent,  which  can  be  evaluated  with  a 
plate  penetration  test)  and  therefore  it  is  advisable  to  treat  ci  and  C2  as  tuning 
parameters,  rather  than  the  sinkage  exponent.  The  influence  of  ci  and  C2  on  slip- 
sinkage  behavior  is  presented  in  Figure  5(b)  where  sinkage  prediction  for  the 
parameters  presented  in  Table  6  are  presented  (the  curves  are  obtained  using  the 
WR  model,  Fz  =  100  N,  and  soil  parameters  from  Table  2). 

3.4.  Shear  Displacement  Calculation 

Another  critical  aspect  of  the  WR  model  is  the  calculation  of  soil  shear  dis¬ 
placement  under  the  wheel.  It  is  beneficial,  at  this  point,  to  define  three  quantities 
that  will  be  discussed  in  this  section: 

•  The  theoretical  tangential  rim  velocity  is  defined  as: 

Vrim  —  fCO  —  r(0(l  —  i)cos0  (13) 

where  r  is  wheel  radius,  co  is  wheel  angular  speed,  i  is  slip,  and  0  is  the 
angular  coordinate  centered  at  the  wheel  axis  (please  refer  to  Figure  11). 
This  is  the  tangential  velocity  of  a  point  on  the  rim  (located  at  angular  po¬ 
sition  0).  It  a  scalar  quantity  and  its  definition  is  derived  from  kinematics. 

•  The  measured  terrain  tangential  velocity  at  the  interface.  This  is  the  ve¬ 
locity  of  the  soil  adjacent  to  the  wheel  rim  in  the  direction  tangential  to 
the  wheel.  It  is  obtained  through  PIV  analysis,  and  is  therefore  a  measured 
quantity  labeled  vpiv.  Being  the  velocity  tangent  to  the  wheel  rim,  this  is  a 
scalar  quantity. 
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•  The  soil  tangential  slip  velocity ,  vt,  which  represents  the  differential  veloc¬ 
ity  between  the  wheel  rim  and  the  soil  along  the  rim  tangential  direction. 
For  measured  data,  vt  can  be  calculated  as  follows: 

Vt  —  vrim  v  piv  (14) 

Wong  and  Reece  postulated  that  soil  shearing  occurs  at  the  wheel-terrain 
interface,  and  that  the  soil  slip  velocity  is  equal  to  the  rim  tangential  velocity 
(for  positive  slip).  Therefore,  they  proposed  the  following  definition  for  the  soil 
tangential  slip  velocity: 


{rco  —  rco(  1  —  i)cosO  +Kvrco(  1  —  i )  [ 6m  <  0  <  6f,i  <  0] 

r(0  —  rco(l  —  i)cos6  [6r  <  0  <  0m.i  <  0]  (15) 

rco  —  rco(l  —  i)cos0  [9r  <  0  <  Of,  i  >  0] 

As  highlighted  in  the  Wheel-Terrain  Interfacial  Stress  section,  the  tangential 
stress  goes  through  a  sign  inversion  for  negative  slip.  This  was  explained  by 
Wong  and  Reece  [24]  as  a  result  of  bulldozed  soil  in  front  of  the  skidding  wheel. 
The  idea  is  that  soil  flows  upward  in  front  of  the  skidding  wheel,  creating  an 
inversion  of  stress.  Hence,  they  modified  the  slip  velocity  for  negative  slip  by 
introducing  a  correction  term  Kv. 

Visual  inspection  of  the  soil  kinematics  under  a  rigid  wheel,  however,  has 
shown  that  soil  at  the  wheel  interface  typically  remains  “attached”  to  the  wheel 
rim,  and  there  is  little  or  no  relative  motion  between  the  soil  and  the  wheel  rim. 
As  a  corollary  to  this,  shearing  does  not  occur  at  the  wheel-terrain  interface,  but 
rather  occurs  along  a  failure  plane  internal  to  the  soil. 

Figure  6  shows  the  theoretical  rim  velocity  together  with  the  measured  soil 
velocity  at  the  rim  (these  are  vectorial  quantities,  not  to  be  confused  with  the 
tangential  soil  velocity  which  is  only  the  component  tangent  to  the  wheel  rim). 
The  data  clearly  show  that  although  the  soil  remains  attached  to  the  wheel,  the 
measured  velocity  deviates  from  the  theoretical  rim  velocity  at  the  entry  and  exit 
angle  where  velocity  gradients  increase. 

Using  PIV  techniques,  it  is  possible  to  quantify,  through  equation  14,  the  slip 
velocity  between  the  wheel  rim  and  the  soil  in  its  proximity.  Figure  7  presents 
results  for  -30%  (a)  and  +30%  (b)  slip.  The  soil  slip  velocity  is  compared  with 
predictions  from  the  WR  model  (equation  15).  For  negative  slip,  the  WR  model 
is  able  to  capture  the  velocity  trend  (see  Figure  7(a)).  The  discontinuity  in  the 
model  is  due  to  the  fact  that  the  slip  velocity  is  defined  as  a  piecewise  equation. 
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(C) 


(d) 


Figure  6:  Soil  velocity  field  for  a  -30%  skidding  wheel  (a)  and  for  a  +30%  slipping  wheel  (b). 
Soil  behavior  is  markedly  different  for  positive  and  negative  slip.  Soil  velocity  at  wheel-soil 
mterface  is  presented  in  figures  (c)  and  (d).  Journal  °f  ^mechanics 
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Figure  7:  Visualization  of  tangential  slip  velocities  for  -30%  (a)  and  +30%  (b)  slip.  vpiV  is  the 
terrain  tangential  velocity  at  the  interface  obtained  through  PIV  measurements.  vt  (Wong  and 
Reece)  is  obtained  through  Equation  15.  vnm  is  the  theoretical  tangential  rim  velocity  presented 
in  Equation  13  and  v,  measured  is  the  soil  tangential  slip  velocity  calculated  through  Equation 
14.  Please  note  that  vt  (Wong  and  Reece)  and  vnm  coincide  for  positive  slip. 

On  the  other  hand,  for  positive  slip  (Figure  7(b)),  the  measured  slip  veloc¬ 
ity  is  not  close  to  the  slip  velocity  predicted  by  WR  model.  (For  positive  slip, 
Wong  and  Reece  postulated  that  the  slip  velocity  corresponds  to  the  rim  veloc¬ 
ity).  This  suggests  that  the  slip  velocity  as  proposed  by  Wong  and  Reece  grossly 
overestimate  the  true  soil  slip  velocity  for  positive  slip. 

However,  it  should  be  noted  that  the  soil  does  not  shear  along  a  thin  shearing 
plane.  In  fact,  a  large  body  of  soil  in  proximity  of  the  rim  (an  area  of  approxi¬ 
mately  10  cm  x  10  cm  in  our  experiments)  is  influenced  by  wheel  passage.  For 
positive  slip,  periodic  formations  of  shearing  bands  are  observed  (approximately 
2  cm  away  from  the  wheel,  as  visible  in  Figure  6(b)),  while  for  negative  slip  it 
is  not  possible  to  identify  clear  shearing  bands.  This  suggests  that  the  Wong  and 
Reece  approach  to  shear  stress  modeling  represents  only  a  coarse  approximation 
to  the  complex  state  of  stress  and  strain  under  the  wheel.  Further  investigation 
into  soil  failure  mechanisms  for  different  slip  levels  is  necessary  to  completely 
understand  how  to  model  shear  stress  under  running  gears. 

3.5.  Terrain  Shear  Modulus  and  Shear  Stress  Modeling 

In  the  WR  model,  the  shear  stress  is  described  as  a  function  of  the  slip  veloc¬ 
ity: 

T  =  (c  +  <Jtan(j))  ( 1  —  e~n  )  (16) 
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j  =  J  vtde  (17) 

where  c  is  the  soil  cohesion,  <7  is  the  normal  stress,  (j)  is  the  angle  of  internal 
friction,  j  is  the  shear  displacement,  and  kx  is  the  soil  shear  modulus. 

This  formulation  is  based  on  two  hypotheses: 

HI  The  soil  slip  velocity  corresponds  to  the  wheel  rim  velocity. 

H2  The  shear  stress  evolution  at  the  interface  is  similar  to  the  shear  vs.  displace¬ 
ment  behavior  observed  during  direct  shear  tests  (or  other  types  of  shear 
tests). 

The  first  hypothesis  has  already  been  discussed  in  section  3.4,  while  the  sec¬ 
ond  hypothesis  will  be  investigated  here. 

Equation  16  was  first  introduced  by  Janosi  and  Hanamoto  [8],  and  was  used 
to  describe  shear  vs.  displacement  behavior  of  granular  materials  in  a  low  den¬ 
sity  state.  The  equation  includes  three  terrain  parameters:  cohesion  c,  angle  of 
internal  friction  (j),  and  shear  modulus  kx.  Cohesion  and  angle  of  internal  fric¬ 
tion  are  well  understood  (and  intrinsic)  soil  parameters,  while  the  shear  modulus 
requires  further  discussion. 

As  we  have  noted  in  [18],  the  true  shear  modulus  kx  (as  presented  in  Table 
2)  for  dry  sand  is  usually  much  smaller  than  what  is  regularly  reported  in  the 
terramechanics  literature.  According  to  the  WR  model,  all  the  shearing  occurs 
within  few  millimeters  from  the  wheel  rim.  Although  it  was  shown  that  this 
is  not  the  case  in  practice,  this  assumption  will  be  held  valid  for  the  sake  of 
discussing  the  WR  model.  Figure  7(b)  suggests  that  the  slip  velocity  for  positive 
wheel  slip  is  highly  over  estimated  using  the  Wong  and  Reece  approach.  As 
a  result,  Figure  8(b)  shows  how  tangential  stress,  predicted  by  the  WR  model, 
is  overestimated  for  positive  slip.  However,  even  for  negative  slip,  where  slip 
velocity  is  in  good  accordance  with  experimental  evidence,  if  the  measured  shear 
modulus  is  utilized  (see  Table  2),  the  resulting  shear  stress  is  overestimated,  as 
shown  in  Figure  8(a). 

The  analysis  of  the  shear  displacement  and  shear  modulus  suggests  that  both 
hypotheses  HI  and  H2  are  not  entirely  true,  and  therefore  the  true  shear  modulus 
kx  is  not  representative  of  soil  shearing  behavior  under  the  wheel.  We  note  that 
shearing  also  depends  on  slip  conditions,  so  kx  is  unlikely  to  remain  constant. 

3.6.  Model  Adaptation 

Thus  far,  the  analysis  has  shown  that  although  classical  terramechanics  mod¬ 
els  may  not  be  able  to  precisely  describe  the  fundamental  phenomena  at  work, 
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(a) 


(b) 


Figure  8:  Comparison  of  Wong  and  Reece  stress  predictions  with  experimental  data  for  -30% 
(a)  and  +30%  (b)  slip.  These  plots  were  obtained  using  a  nominal  kx  =  0.0006  and  other  soil 
parameters  presented  in  Table  2.  Tangential  stress  is  overestimated  for  both  positive  and  negative 
slip.  This  is  because  the  nominal  kx.  Note  that  for  positive  slip,  there  is  a  shift  between  WR 
model  and  measurements:  this  is  because  for  that  specific  combination  of  slip  and  normal  load, 
WR  model  is  underestimating  sinkage  (i.e.,  entry  angle  is  smaller).  Exhaustive  analysis  of  WR 
model  behavior  is  presented  in  Figures  9,  10,  and  Table  7. 


they  are  still  able  to  describe  all  the  salient  characteristic  of  wheel  terrain  interac¬ 
tion.  Unfortunately,  these  methods  remain  empirical  in  nature,  and  require  some 
tuning  in  order  to  deliver  satisfactory  results.  However,  the  analysis  presented  in 
this  paper  has  led  us  to  conclude  that  it  is  not  necessary  to  introduce  additional 
coefficients  to  produce  reasonable  model  predictions. 

Considering  the  results  obtained  from  the  analysis  of  shear  displacement  and 
shear  modulus,  we  propose  to  retain  the  Wong  and  Reece  shear  displacement 
model  untouched  and  modify  the  definition  of  shear  modulus  kx.  In  fact,  if  the 
shear  modulus  kx  is  made  function  of  slip,  it  is  possible  to  accurately  model 
wheel  behavior  for  a  relatively  large  range  of  vertical  loads,  and  for  slip  up  to 
±70%.  (  We  limit  our  analysis  to  ±70%  since  our  experimental  data  was  col¬ 
lected  within  this  range).  We  argue  that  the  shear  modulus  should  be  viewed 
as  an  empirical  parameter  that  governs  the  complex  mapping  between  soil  flow 
and  shear  stress  at  the  wheel  interface,  and  this  study  has  shown  that  it  cannot 
be  simply  extrapolated  from  direct  shear  tests  and  assumed  a  constant  parameter 
value. 

Figure  9  shows  a  comparison  of  torque  M,  traction  Fx,  and  sinkage  z  between 
experimental  data  and  the  WR  model.  For  the  WR  model,  the  parameters  con¬ 
tained  in  Table  2  and  Table  6  (the  first  two  columns)  have  been  used,  with  the 
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exception  of  kx,  which  has  been  modified  as  follows: 

kx  =  cx{i)kx  o  (18) 

where  cx  is  a  correction  factor  (a  function  of  wheel  slip)  and  kx o  is  the  nominal 
shear  modulus  (as  presented  in  Table  2). 

The  correction  factor  was  found  by  minimizing  the  following  cost  function: 

f  =  (wr2  zde-MSA^j  (19) 

where  MSA  is  the  torque  measured  through  the  sensing  array  and  the  integral  term 
is  the  torque  as  predicted  by  the  Wong  and  Reece  model  (see  APPENDIX  A  for 
details).  Parameters  for  determining  the  correction  factor  cx  can  be  found  in 
Table  8. 

The  piece-wise  linear  trend  of  the  correction  factor  cx  can  be  explained  in 
light  of  the  failure  mechanism  under  the  wheel. 

For  large  negative  slip,  braking  force  is  primarily  induced  by  normal  stress, 
while  shear  stress  has  relatively  less  influence.  In  fact,  Figure  3  shows  that  shear 
stress  is  approximately  three  times  smaller  for  -70%  slip  as  compared  to  +70% 
slip.  Because  of  this,  the  shear  modulus  grows  for  negative  slip,  as  presented  in 
Figure  9(d).  On  the  other  hand,  for  small  positive  slip,  the  soil  under  the  wheel 
does  not  develop  full  shearing  bands,  and  therefore  the  shear  modulus  must  be 
large  in  order  to  model  the  lack  of  traction  produced.  For  increasing,  positive 
slip,  clearly  discernible  shear  bands  develop  in  the  soil,  and  therefore  the  shear 
modulus  again  approaches  its  nominal  value  (i.e.,  the  correction  factor  goes  to 
one). 

An  alternative  approach,  aimed  at  simplifying  the  model,  is  to  utilize  the 
same  shear  velocity  formulation  for  both  positive  and  negative  slip: 

vt  —  rco  —  rco(l  —  i)cos0  (20) 

As  shown  in  Figure  7  this  is  not  strictly  correct,  however,  considering  the 
empirical  nature  of  this  approach,  it  is  not  unreasonable.  Note  that  this  simplifi¬ 
cation  only  affects  negative  slip  while  the  model  remains  unchanged  for  positive 
slip.  The  main  benefit  of  this  approach  is  that  the  model  becomes  continuous 
across  the  zero  slip  abscissa. 

The  simplified  model  is  compared  to  experimental  data  in  Figure  10.  Pre¬ 
diction  of  torque,  drawbar  force,  and  sinkage  with  the  simplified  shear  model  is 
close  to  what  was  obtained  with  the  original  model  (it  is  identical  for  positive 
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Figure  9:  Comparison  of  Wong  and  Reece  model  with  experimental  data  (sensing  array)  utilizing 
a  correction  factor  (d)  for  the  shear  modulus  kx. 


26 


Journal  of  Terramechanics 


30 


C.  Senatore  and  K.  Iagnemma 


(a) 


(b) 


(c)  (d) 

Figure  10:  Comparison  of  the  simplified  model  with  experimental  data  utilizing  a  correction 
factor  for  the  shear  modulus  kx.  The  simplified  model  produces  identical  results  to  Wong  and 
Reece  model  for  positive  slip  and  overall  similar  results  for  negative  slip.  The  main  advantage  of 
the  simplified  model  is  to  simplify  the  definition  of  shear  displacement,  avoiding  discontinuity 
across  zero  slip. 
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Table  7:  Absolute  error  A  and  full  scale  percent  error  E/  of  the  WR  model  with  respect  to  the 
sensing  array  readings.  Average  error  stays  below  11%  with  more  than  half  of  the  data  points 
below  5%  error. 
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slip).  It  is  interesting  to  note  that  the  shear  modulus  correction  factor  assumes  an 
exponential  decay.  The  reasons  behind  this  behavior  are  the  same  as  those  high¬ 
lighted  above  for  the  original  model.  However,  with  the  simplified  model,  the 
correction  factor  becomes  continuous  because  the  shear  displacement  equation 
is  continuous.  The  highly  nonlinear  behavior  of  cx  for  negative  slip  is  caused  by 
the  fact  that  with  the  simplified  model,  the  shear  stress  does  not  undergo  a  sign 
inversion.  This  causes  the  torque  to  grow  larger  for  negative  slip,  and  therefore 
it  requires  a  larger  cx,  while  in  the  original  model  the  integral  of  r  was  internally 
reduced  by  the  sign  change. 

As  shown  in  Table  7,  the  model  is,  on  average,  accurate  within  11%.  For  a 
balanced  analysis  it  is  important  at  examine  the  full  scale  error  £/-.  Sinkage  error, 
although  large  in  a  relative  sense  is  typically  on  the  order  of  2-5  millimeters.  This 
is  unsurprising  considering  the  practical  difficulty  of  preparing  a  perfectly  flat 
soil  surface  for  testing.  In  practice,  run-to-run  variation  on  the  order  of  several 
mm  is  impossible  to  avoid. 

Table  8:  Correction  factor  parameters  for  the  Wong  and  Reece  model  and  for  the  simplified 
model. 


Wong  and  Reece  (cx  —  pii  +  pi) 

Simplified  ( cx  —  p\eP21) 

i  <  0 

i  >  0 

- 

Pi 

-31 

-14 

20.4 

P2 

-3 

16 

-2.9 

4.  Conclusions 

This  paper  presented  an  analysis  of  the  interaction  phenomena  governing 
lightly  loaded  rigid  wheel  performance  on  dry  sand.  The  stress  analysis  high¬ 
lighted  that  stress  distribution  across  wheel  width  varies  significantly  only  for 
large  slip.  For  large  positive  slip,  stress  is  higher  at  the  wheel  edges,  while 
for  large  negative  slip  the  opposite  is  true.  Comparison  with  the  Wong  and 
Reece  model  showed  that  this  model  is  in  theory  able  to  characterize  the  mo¬ 
bility  performance  of  lightweight  vehicles.  However,  the  empirical  nature  of 
method  should  not  be  forgotten.  It  was  found  that  the  coefficients  c\  and  C2  do 
not  depend  on  vertical  load,  and  therefore  can  be  considered  constant  for  a  spe¬ 
cific  wheel-soil  configuration.  The  analysis  of  soil  kinematics  under  the  wheel 
showed  that  the  hypotheses  behind  the  shear  stress  formulation  are  not  entirely 
valid:  this  explains  why  the  shear  modulus  kx,  obtained  from  direct  shear  tests 
(or  other  shear  tests),  does  not  produce  accurate  results  when  used  in  the  Wong 
and  Reece  model. 
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APPENDIX  A 

The  Bekker-Wong-Reece  model  relies  on  analysis  of  two  fundamental  rela¬ 
tionships:  the  pressure-sinkage  relationship,  and  the  shear  stress-shear  deforma¬ 
tion  relationship.  In  the  context  of  rover  mobility,  the  pressure-sinkage  relation¬ 
ship  governs  the  depth  that  a  rover  wheel  will  sink  into  the  terrain  and  therefore 
how  much  resistance  it  will  face  during  driving.  The  shear  stress-shear  displace¬ 
ment  relationship  governs  the  amount  of  traction  that  a  wheel  will  generate  when 
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driven  and  therefore  how  easily  it  will  progress  through  terrain  and  surmount  ob¬ 
stacles.  The  pressure- sinkage  relationship  was  described  by  Bekker  in  the  form 
of  a  semi-empirical  equation  that  relates  sinkage  with  normal  pressure  of  a  plate 
pushed  down  into  the  soil.  The  proposed  relation  is  commonly  referred  as  the 
Bekker  equation,  and  provides  a  link  between  the  kinematics  (sinkage)  and  stress 
(pressure)  of  a  plate  (which  can  be  viewed  as  a  proxy  for  a  wheel  or  track): 


(21) 


Parameters  kc,  kQ,  n  are  empirical  constants  that  are  dependent  on  soil  prop¬ 
erties,  while  b  corresponds  to  the  plate  width.  These  parameters  can  be  obtained 
from  field  tests  conducted  with  a  device  called  a  bevameter.  Specifically,  stress 
can  be  divided  in  two  components  (assuming  a  two  dimensional  model,  and  mo¬ 
mentarily  ignoring  out  of  plane  motion):  normal  stress  and  tangential  stress.  A 
schematic  representation  of  the  stress  distribution  at  a  wheel-terrain  interface  is 
presented  in  Figure  1 1 .  Normal  stress  can  be  calculated  by  starting  with  Bekker’s 
pressure- sinkage  relation,  and  introducing  a  scaling  function  intended  to  satisfy 
the  zero-stress  boundary  conditions  present  at  the  fore  and  aft  points  of  contact  of 
the  wheel  with  the  terrain  (known  as  “soil  entry”  and  “soil  exit”).  The  equation 
is  expressed  as  a  piece- wise  function,  as  follows: 


a 


Figure  1 1 :  Schematic  representation  of  normal  and  tangetial  stress  profile  along  a  the  wheel-soil 
interface. 
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j 0i  =  (f  +  k<i>) zi  e»'  <°<of 
y  02  —  (jfr  +  j  Or  <  6  <  Om 
Zi  —  r(cosO  —  cos  Of) 

Z2  =  r  ^cos  ^/--^—^-(Of-O,,,)^ -cosO/^j  (22) 

where  Of  is  the  soil  entry  angle,  0,  is  the  exit  angle,  and  0m  is  the  angle  at  which 
the  maximum  normal  stress  occurs  (see  Figure  11).  Wong  and  Reece  suggested 
use  of  a  linear  equation  to  calculate  0m : 

9m  =  (C|  +C2l)0f  (23) 

The  shear  stress-shear  displacement  relationship  is  based  on  the  Mohr-Coulomb 
failure  criterion,  coupled  with  a  modulation  function  proposed  by  Janosi  and 
Hanamoto  [8]: 


z  —  (c  +  otan(j))  —  e  ^  (24) 

where  c  is  the  soil  cohesion,  (j)  is  the  angle  of  internal  friction,  kx  is  the  shear 
modulus  (a  measure  of  shear  stiffness),  and  j  is  shear  deformation: 

rr o  ref  dO 

j=  /  vtdO=  /  vr—  (25) 

Jo  Je  (o 

where  vt  is  the  tangential  slip  velocity  and  kx  is  the  shear  modulus.  Expressing 
vt  as  a  function  of  angular  velocity  and  slip,  it  is  possible  to  obtain  an  expression 
for  shear  displacement  in  closed  form  for  positive  slip: 

vt  —  rco  —  r(o{\  —  i)cosO  (26) 

j+  =  r  [Of  —  0  —  (1  —i)(sinOf  —  sinO)]  (27) 

For  negative  slip,  starting  from  the  observation  that  soil  flows  upward  at  the 
leading  edge  of  the  wheel,  Wong  and  Reece  proposed  a  modified  equation  for 
the  slip  velocity  between  Of  and  0m: 

vt  =  r(0  —  r(o{\  —  i)cosO  +  Kvr(o{\  —i)  0m  <  0  <  Of  (28) 

The  constant  Kv  was  calculated  by  requiring  that  the  shear  displacement 
should  be  zero  for  0  —  Of  —  0m,  leading  to: 
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Kv  = 


1 

(i-0 


(1  —  /)  (sin  Of  —  sin  0m) 
Of  —  Om 


(29) 


It  should  be  noted  that  Wong  and  Reece  proposed  a  different  approach  for 
calculating  0m  for  negative  slip,  while  here  we  will  assume  that  the  same  linear 
expression  can  be  used.  Also,  Wong  and  Reece  used  positive  values  for  skid; 
hence  the  sign  of  slip  i  is  inverted  here  (slip  is  defined  by  Equation  5).  Therefore, 
the  shear  displacement  for  negative  slip  can  be  written  as: 


j 


r[(0/-0)(l+^v(l-/))-(l-/)(sin0/-sin0)]  8m  <  0  <  0/ 
r(0m  —  0  —  (1  —  z)(sin 0m  —  sin 0))  8r  <  0  <  Bm 


This  definition  of  shear  displacement  allows  the  Wong  and  Reece  model  to 
predict  shear  displacement  sign  inversion.  Once  the  stress  profile  acting  on  a 
wheel  has  been  completely  defined,  these  profiles  can  be  integrated  to  determine 
the  net  forces  and  torques,  which  are  then  summed  to  compute  overall  vehicle 
motion.  Traction  forces  generated  by  a  wheel  can  be  decomposed  in  two  com¬ 
ponents:  a  thrust  component,  which  acts  to  move  the  vehicle  forward;  and  a 
compaction  resistance  component,  which  resists  forward  motion.  Thrust,  T,  is 
computed  as  the  sum  of  all  shear  force  components  in  the  direction  of  forward 
motion: 

f0f 

T  —  wr  tcos  Odd  (31) 


Compaction  resistance,  Rc,  is  the  result  of  all  normal  force  components  act¬ 
ing  to  resist  forward  motion,  and  can  be  thought  of  as  the  net  resistance  force 
provided  by  the  soil: 

rdf 


Rc  —  wr  (7  sin  Od  0 


(32) 


The  net  longitudinal  force,  also  termed  the  drawbar  pull,  Fx  ,  is  calculated  as 
the  difference  between  the  thrust  force  and  resistance  force.  Fx  is  the  resultant 
force  that  can  provide  a  pulling/braking  force  at  the  vehicle  axle. 


Fx  —  T  —  Rc 


(33) 


The  importance  of  drawbar  force  is  obvious,  since  a  positive  drawbar  force 
implies  that  a  rover  can  generate  forward  motion  on  a  particular  patch  of  terrain, 
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while  a  negative  drawbar  force  suggests  that  forward  motion  is  difficult  or  im¬ 
possible.  Torque,  M,  is  the  resultant  of  shearing  action  along  wheel  rim,  and  can 
be  calculated  as: 


M  =  wr 


T  dO 


(34) 


The  sinkage  of  a  wheel  can  be  calculated  by  solving  a  vertical  force  equilib¬ 
rium  problem,  which  enforces  the  fact  that  the  force  resisting  wheel  penetration 
into  the  soil  must  be  balanced  by  the  vertical  load  acting  on  that  wheel. 


Fz  —  wrf  rsin  ddB  +  wrf  a  cosOd  6  (35) 

J  0r  J  6r 

Please  note  that  the  balance  of  vertical  load  is  the  first  equation  that  must  be 
solved  in  order  to  define  entry  and  exit  angles  Of  and  0r  (as  previously  noted, 
exit  angle  0 ,  is  usually  held  constant). 
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